Take-home Exercise 3: NKDE

Packages

pacman::p_load(sf,tidyverse,spNetwork,tmap,classInt, viridis)

Import data

crime_data <- read_rds("../../data/proj_data/rds/accidents_thai.rds") %>% 
  st_transform(crs = 24047)
# thai_bound <- st_read(dsn = "../../data/proj_data/thai_adm_boundary",
#                       layer = "tha_admbnda_adm1_rtsd_20220121") %>% 
#   st_transform(crs = 24047)
thai_bound <- st_read(dsn = "../../data/proj_data/thai_adm_boundary",
                      layer = "tha_admbnda_adm2_rtsd_20220121") %>%
  st_transform(crs = 24047)

Reading layer `tha_admbnda_adm2_rtsd_20220121' from data source 
  `/Users/tangtang/Desktop/IS415 Geospatial Analytics and Applications/practice/is415gaa/data/proj_data/thai_adm_boundary' 
  using driver `ESRI Shapefile'
Simple feature collection with 928 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 97.34336 ymin: 5.613038 xmax: 105.637 ymax: 20.46507
Geodetic CRS:  WGS 84

thai_road <- read_rds("../../data/proj_data/rds/thai_roads.rds") %>% 
  st_transform(crs = 24047)

Visualising difference in network constrained kernel density through different filters by columns

Some potential filtering to apply to dataset to allow user to zoom in into specific area he/she is interested in:

• Accident categories

• Fatal accident [yes/no]

• Province (MUST-HAVE)

But let’s visualize the overall look

plot(crime_data)

Example: filter by province = “Bangkok”, accident_categories = “all”, fatal_accident == “all” (Allowing users to change freely)

Narrow down to city due to computation limitation

province_factor = "Bangkok"

crime_data_filtered = crime_data %>% 
  filter(province_en == province_factor)

thai_prov <- thai_bound %>% 
  filter(ADM1_EN == province_factor)

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(thai_prov) +
  tm_polygons() +
tm_shape(crime_data_filtered) +
  tm_dots(col = "black",
          size = 0.01)

tmap_mode("plot")

tmap mode set to plotting

Network-constrained Kernel Density Estimation

Get the portion of road network for the selected town from selected province

thai_road <- thai_road %>%
  filter(fclass %in% c("motorway","trunk","primary","secondary","tertiary"))
road_prov <- st_intersection(thai_road, thai_prov)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

road_prov_sf <- st_cast(road_prov,"LINESTRING")

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

Warning in st_cast.MULTILINESTRING(X[[i]], ...): keeping first linestring only

any(duplicated(road_prov_sf))

[1] FALSE

Prepare lixel object

1st try with 1000m diameter

lixels <- lixelize_lines(road_prov_sf,lx_length = 1000,mindist = 500)

Generate line centre points

Likely to allow adjustment of kernel type, bandwidth set to be equivalent to mindist of lixels

samples <- lines_center(lixels)

densities <- nkde(road_prov_sf,
                  events = crime_data_filtered,
                  w = rep(1, nrow(crime_data_filtered)),
                  samples = samples,
                  ## important
                  kernel_name = "quartic",
                  bw = 500,
                  ## important
                  div = "bw",
                  method = "simple",
                  digits = 1,
                  tol = 1,
                  grid_shape = c(1,1),
                  max_depth = 8,
                  # aggregate events within 10m radius (faster calculation)
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE
                  )

Generate plot

samples$density <- densities
lixels$density <- densities

samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

tmap_mode('view')

tmap mode set to interactive viewing

tm_shape(lixels) +
  tm_lines(col="density",
           lwd = 3)

tmap_mode('plot')

tmap mode set to plotting

UI design for Analysis

Storyboard of NKDE

The user would be able to filter the dataset to see impact different categories on the network-constrained kernel density map created and use the interpretation rules given to derive at a conclusion.